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We have developed a numerically exact approach to 
compute real-time path integral expressions for quantum 
transport problems out of equilibrium. The scheme is 
based on a deterministic iterative summation of the path 
integral (ISPI) for the generating function of nonequi- 
librium observables of interest, e.g., the charge current 
or dynamical quantities of the central part. Self-energies 
due to the leads, being nonlocal in time, are fully taken 
into account within a finite memory time, thereby in- 
cluding non-Markovian effects. Numerical results are 
extrapolated first to vanishing (Trotter) time discretiza- 
tion and, second, to infinite memory time. The method 
is applied to nonequilibrium transport through a single- 
impurity Anderson dot in the first place. We bench- 
mark our results in various regimes of the rich param- 
eter space. In the respective regime of validity, ISPI 
results are shown to match those of other state-of-the 
art methods. Especially, we have chosen the mixed va- 
lence regime of the Anderson model to compare ISPI 
to time dependent density matrix renormalization group 
(tDMRG) and functional RG calculations. Secondly, we 
determine the nonequilibrium current I{V) through a 
molecular junction in presence of a vibrational mode. 
We have found an exact mapping of the single impurity 
Anderson-Holstein 



model to an effective spin-1 problem. In analytically 
tractable regimes, as the adiabatic phonon or weak 
molecule-lead coupling regime, we reproduce known 
perturbative results. Studying the crossover regime be- 
tween those limits shows that the Franck-Condon block- 
ade persists in the quantum limit. At low temperature, 
the Franck-Condon steps in the I{V) characteristics are 
smeared due to nonequilibrium conditions. The third 
system under investigation here is the magnetic An- 
derson model which consists of a spinful single-orbital 
quantum dot with an incorporated quantum mechani- 
cal spin- 1/2 magnetic impurity. Coulomb interaction to- 
gether with the exchange coupling of the magnetic impu- 
rity with the electron spins strongly influence the dynam- 
ics. We investigate the nonequilibrium tunneling cur- 
rent through the system as a function of exchange and 
Coulomb interaction as well as the real-time impurity 
polarization. From the real-time evolution of physical 
observables, we are able to determine characteristics of 
the time-dependent nonequilibrium current and the re- 
laxation dynamics of the impurity. These examples illus- 
trate that the ISPI technique is particularly well suited 
for the deep quantum regime, when all time and energy 
scales are of the same order of magnitude. 
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1 Introduction Small condensed matter systems that 
behave as transistors have attracted numerous research ac- 
tivities in recent years. When placed in a transport setup, 



current and noise measurements are used to determine the 
properties of the usually strongly correlated quantum me- 
chanical system. In as well as out of thermal equilibrium, 
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several quantum many-body properties of systems such as 
quantum dots or single molecule junctions are accessible 
in experiments |1,2|. By down-scaling electronic transis- 
tors to the nanometer scale, quantum dots (QD's) allow 
to control electronic and/or spin properties of single elec- 
trons and might be used to impress and read out informa- 
tion. Single molecules, inheriting also vibrational degrees 
of freedom, are ideal candidates for such devices since 
they can rather easily be functionalized such that current 
switches are implemented in a nanoscale environment fT 
[4 , 5 6 7 8 1 . Furthermore, a controlled design of functional 
groups attached to molecules seems in reach. The field 
of molecular electronics is at the intersection of interdis- 
ciplinary research areas, combining chemical and physi- 
cal properties of molecules. Under nonequilibrium condi- 
tions transport through molecules challenges experimen- 
talists as well as theoreticians. Nonequilibrium in this re- 
view is referred to as the situation that the Fermi ener- 
gies, provided by source and drain electrodes and located 
'left' and 'right' of the structure, are not equal such that 
electrons are transferred through the system by tunneling. 
The voltage is assumed to drop at the nanostructure. In 
such finite voltage bias situations, many interesting physi- 
cal effects arise due to the quantum nature of the electrons 
and their strong Coulomb interaction among each other 
and/or by their interaction with molecular vibrations or lo- 
cal magnetic moments. Prominent examples are resonant 
tunneling or the Kondo effect ||91|10II11L Different tech- 
niques are used to approach the strong coupling limit. For 
instance, the regime of low energy, low temperature/bias 
has been approached by Fermi liquid theory |12|, via in- 
terpolative schemes |13], using integrability concepts IT4l . 
or by the perturbative renormalization group (RG) jfTSl 
[T6l . A perturbative RG analysis has been performed by 
Schoeller and Konig 1 17 1. The nonequilibrium generaliza- 
tion of Wilson's numerical RG approach lITSi and of the 
fRG method |19| have been discussed. Transport features 
have also been discussed by perturbation theory in the in- 
teraction strength |20|. On the other hand, perturbatively 
treating the tunneling matrix elements is a powerful way of 
describing the incoherent regime, see Ref. [9|. Density ma- 
trix renormalization group techniques have been extended 
to the nonequilibrium regime |21 22|. In addition the flow 
equation method allows to study the Toulouse point of the 
Kondo model |23|. 

In case that a molecule is placed between electrodes, an 
appropriate theoretical description is needed to deal with 
its characteristic vibrational (phonon) degrees of freedom 
Il2l l24i . The interplay between mechanical and electronic 
degrees of freedom is of interest in other areas of physics 
as well, e.g., for inelastic tunneling spectroscopy ||251 . na- 
noelectromechanical systems |26|, break junctions ll27l . 
and suspended semiconductor or carbon-based nanostruc- 
tures ||281|29II301I3TI . Including vibrational degrees of free- 
dom via a simple Anderson Holstein (AH) model, where 
a spinless electronic level is coupled to a single oscillator 



mode, shows several effects as the Franck-Condon block- 
ade, negative differential conductance, or current induced 
heating or cooling |32,33|; for a review, see Ref. |25|. As 
for the interacting Anderson model, analytical approaches 
typically address different comers of parameter space, a 
full theory that connects those corners seems not in reach 
at present. 

A third topic adressed in this review is the investiga- 
tion of a magnetic QD. Those setups have been studied 
experimentally in ensembles which are particularly suited 
for the investigation by laser and electromagnetic fields 
Il34||35|[36l[37i r38 39 40 1 . They are designed with standard 
lithographic methods and are technologically well estab- 
lished. Moreover, embedding individual magnetic Mn ions 
into quantum dots and studying the electrical properties is 
possible II41II42||43| . Small quantum dots with few charge 
carriers and a single magnetic impurity may become im- 
portant candidates for efficient high density spintronic de- 
vices. 

There is a considerable need for numerical methods 
which describe small quantum systems out of equilibrium 
accurately and, ideally, treat simple model Hamiltonians 
numerically exactly. Numerical renormalization group ll44ll 
or quantum Monte Carlo (QMC) calculations l!45lH6ll47l 
48 , 49 , 50 1 provide a possible line of attack to those prob- 
lems. Due to the dynamical sign problem, these calcula- 
tions become increasingly difficult at low temperatures, 
but in several parameter regions, the stationary steady-state 
regime seems accessible. Based on non-standard ensem- 
bles, the steady state is described in a recent work by Han 
lISTI . using an imaginary-time QMC approach followed by 
a double analytical continuation scheme. This last step is 
numerically the most difficult part |52|. 

We here review the development of a novel numeri- 
cal scheme denoted as iterative summation of real-time 
path integrals (ISPI), in order to address quantum trans- 
port problems out of equilibrium B53l . Many-body systems 
driven out of equilibrium are known If l5ll54ll55l to acquire 
a steady state that may be quite different in character from 
their ground-state properties. Details of the steady state 
may depend on the nature of the correlations, as well as on 
the way in which the system is driven out of equilibrium. 
Our ISPI approach, described in detail below, provides an 
alternative and numerically exact method to tackle out-of- 
equilibrium transport in correlated quantum dots. Based on 
the evaluation of the full nonequilibrium Keldysh generat- 
ing function, along with the inclusion of suitable source 
terms, observables of interest are computed. It builds on 
the fact that nonlocal in time correlations, induced by the 
fermionic leads, decay exponentially in the long-time limit 
at any finite temperature. Within a characteristic time Tc, 
all correlations are taken into account, while for larger 
times, the correlations are dropped due to their exponen- 
tially small contributions. This allows us to construct an 
iterative scheme to evaluate the generating function. An 
appropriate extrapolation procedure allows to eliminate the 
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Trotter time discretization error (the Hubbard- Stratonovich 
(HS) transformation below requires to discretize time) as 
well as the finite memory-time error. This yields the de- 
sired numerically exact value for the observables of inter- 
est. Note that the need for a finite memory time makes our 
approach difficult to apply at very low energies (T, V — > 
0). Fortunately, other methods are available in this regime. 
At finite T or V, the requirement of not too long memory 
times is exploited and the spin path summation remains 
tractable. Recently, also Segal et al., see Ref. 1*561 have 
provided an alternative formulation of the ISPI approach 
in terms of Feynman Vernon like influence functionals. 

The ISPI scheme is implemented here for three char- 
acteristic impurity models. First, the single-level Ander- 
son impurity model |57,58,59 60 61 1, second the spinless 
Anderson Holstein model |25|, which mimics the behav- 
ior of a molecular quantum dot. Finally, when magnetic 
molecules are investigated, additional couplings to local- 
ized spin impurities are important. Results for the relax- 
ation dynamics of such a system are presented in this arti- 
cle as well. 

The present paper is organized as follows. In Sec. |2] 
we introduce the model for the quantum dot coupled to 
normal leads. We present the computation of the generat- 
ing function for the nonequilibrium Anderson model there 
as well. The presence of an external source term allows to 
calculate the current as a functional derivative. In Sec. [3] 
we introduce the numerical iterative path integral summa- 
tion method, from which we obtain observables of interest. 
We give a detailed discussion of the convergence proper- 
ties of our method and describe the extrapolation scheme. 
For several sets of parameters, we will present results and 
benchmark checks in Sec. |4] In addition, results for the 
mixed valence regime of the Anderson model are pre- 
sented. A good agreement with tDMRG and fRG is found. 
The Anderson Holstein model is studied by means of the 
ISPI method in Sec. E\ the main finding of a sustained 
Franck-Condon blockade at low temperatures is discussed 
upon validating our spin-1 mapping of the AH Hamilto- 
nian. By far the most complex model of this work, a mag- 
netic interacting QD is studied in Sec|6] We briefly review 
the necessary changes in the appearing Keldysh generating 
function and discuss our results on the relaxation time and 
impurity polarization as well as the tunneling current. A 
summary is given in Sec. U\ 

2 Keldysh generating function for non-interacting 
impurity models We consider the generic Hamiltonian 
for a quantum dot that is coupled to metallic leads to the 
left and right side (h — l,kB = I) 

"H = Hdot + Hieads + Ht 



kpa 



/ . Y'P^kp(T^<^ 



(1) 



Here, Eq^ = _Eo + <jB with a =t, \-= ± is the energy of 
a single electron with spin a on the isolated dot. Tuning a 
back gate voltage or a Zeeman magnetic field term oc B 
changes the value of E'ocr. The latter is assumed not to af- 
fect the electron dispersion in the leads. The corresponding 
dot electron annihilation/creation operator is da/d\, with 
the density operator n^ = d\da- Possible eigenvalues of 
Ug. are J^ = 0, 1, corresponding to the empty or occupied 
electronic state with spin a. Interactions on the quantum 
dot are treated in Sec. 12. II and not considered for the mo- 
ment. In Eq. (fTb, efep denotes the energies of the noninter- 
acting electrons (operators Ckpa) in lead p = L/R = ±, 
with chemical potential fip ~ peV/2. Quantum dot and 
leads are connected by the tunnel couplings tp. The ob- 
servable of interest is the (symmetrized) tunneling current 



m 



E hp( 



j^kpa/t 



J 



-ptp{(''kpja)t 



kpa 



(2) 



where Ip{t) = -eNp{t) with Np{t) = {Eka4pa(^kpa)t- 
The stationary steady-state dc current follows as the 
asymptotic long-time limit, / = limi_j.oo lit)- We have ex- 
plicitly confirmed that current conservation, I^ + Ir — 0, 
is numerically fulfilled for the ISPI scheme. 

In the presence of a finite bias voltage, V ^ {), the 
Keldysh technique 1 62, 63, 64 1 provides a way to study 
nonequilibrium transport. In this formalism, the time 
axis is extended to a contour with a — ± branches, see 
Ref. 1 64 1, along with an effective doubling of fields. The 
Keldysh Green function (GF) exhibits a matrix structure 

Gtfito.,t'i^) = -i{Tc[Mt^)^}{t'p)]), where Tc denotes 
the contour ordering of times along the Keldysh contour, 
and i, j — L,R,0 correspond to fields representing lead or 
dot fermions, respectively. We omit the spin indices here, 
remembering that each entry still is a diagonal 2x2 matrix 
in spin space. The Keldysh partition function contains all 
relevant information about the physics of the system. In 
order to obtain it, we first integrate over the noninteracting 
lead fermion fields. Subsequently, we integrate over the 
dot fields as well. In a fermion coherent state basis, the 
generating function is 



zM = / 2? 



]^da,d<T, 



Ckpa 1 ^kpa 



^iS[d(j ,da ,Ckp^ ,Ckp^] 



(3) 

with Grassmann fields {da-,da,c,c). The external source 
term, which allows to compute the current at measurement 
time tm, is chosen such that 



/(t. 



'*^lnZ[ry] 



(4) 



r(=0 



kpa 



Correspondingly, it is also possible to evaluate other ob- 
servables, e.g., the zero-frequency shot noise, by introduc- 
ing appropriate source terms and performing the corre- 
sponding derivatives. The action is S" = Sdot + Sieads + 
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St + Srj, see Ref. Il53l for the explicit expressions. After 
integrating over the leads' degrees of freedom, the effective 
action for the dot becomes nonlocal in time. The generat- 
ing function for the noninteracting system reads 



with 



JJdcrda 



^i{Sdot,0 + Senv) 



(5) 



Senv^ I dt I dt'^d,{t)[-iL{t,t')+-fR{t,t') 

+ ^[7L(t,t')-7i..(i,0] 

X [5{t - tra) + 5{t' - tra)]\d,{t'). 



(6) 



For the source term, the physical measurement time t„i is 

fixed on the upper (+) branch, hence the ( ) Keldysh 

element of the source term self-energy vanishes. The 
7p(t, i') matrices in Eq. ^ represent the leads, their 
Fourier transforms in frequency space are explicitly given 
as 2 X 2 Keldysh matrices 

With the usual assumption that the leads are in thermal 
equilibrium, /(w) = \/{e'^/"^ + 1). Taking the wide- 
band limit with a constant density of states p{ep) per 
spin channel around the Fermi energy, the hybridization 
Fp = 7rp(ei?)|ipp of the dot level with lead p enters. We 
focus on symmetric contacts Fj^ = Fji = F/2 and on sym- 
metrically applied bias voltages as well. The generalization 
to asymmetric contacts is straightforward. 

In the next step (still for vanishing on-dot interactions), 
we integrate over the dot degrees of freedom. This yields 
the noninteracting generating function 

Zm[v\ = n^ct [lG^^(t,t')+T^S\t,t')\ . (8) 
a 

The function G^lit, t') follows from 

Goct(w) = [{uj - eo„)T:, - 7l(w) - 7_R(a;)]"^ 
1 



1 + [(c^ - eo,)/F]^ 

fio-eoa + iF{l-F) 
[ iF{F-2) 



iFF 

-Lu + eoa + iFil-F)^ 

where t^ is the standard Pauli matrix in Keldysh space, 
and F = f{uj + eV/2) + f{uj - eV/2). Moreover, the 
self-energy for the source term is obtained as 

X [S{t - tra) + Sit' - t,„)] . (9) 

Up to this point, we have discussed the noninteracting case. 
The next section describes how to include the electron- 
electron interaction. 



2.1 Hubbard-Stratonovich transformation In the 

presence of on-dot Coulomb interactions, we add the 
Coulomb term 

Hint = Un-fHi (10) 

to the Hamiltonian in Eq. ([TJ. For our purpose, it is conve- 
nient to use the operator identity n^ni ~ ^{n^ + n^) — 
i(n^ — n^)^, which results in a shift of the single-particle 
energies eoa = Eoa + U/2, and we may rewrite Hdot = 
Hdotfl + Hu = Y.a ''O^na - f (n^ - n;)^. Now, the ac- 
tion in Eq. p} in real time contains quartic terms of dot 
Grassmann fields and a Gaussian integration is not pos- 
sible. We will use a time-discrete path integral for the 
following discussion ll65l . In order to decouple the quar- 
tic term, we discretize the full time interval, t = NSt, 
with the time increment 6t- On each time slice, we per- 
form a Trotter breakup of the dot propagator according to 

piSt(Ho+HT) _ ^IStHr /2 ^iStHo giSfHT /2 _|_ QiS'^) where 

Hq = Hdot + Hieads- According to Refs. ||g6l|g71l68l the 
emerging Trotter error can be systematically eliminated 
from the results |69,70|, see below. On a single Trotter 
shce, a discrete Hubbard-Stratonovich transformation ||67l 
I7ni72ll73l allows to decouple the interaction. This locally 
introduces Ising-like discrete spin fields s„ — (s+, s^) on 
the a = ± branches of the Keldysh contour with s" — ±1 
on the n-th Trotter slice. For a given Trotter slice, we define 



^±iStU{n-f-7i^f/2 



1 y^ g-5tA±s±(n^-«^)^ 
s±=± 



(11) 



The HS parameter is obtained from the equation, 

cosh((5tA±) = cos{6tU/2) ±isin{6tU/2), 



under the condition that U > 0, see Ref. 1531 . Note 
that the arbitrarily chosen overall sign of A± does not 
influence the physical result. Uniqueness of this HS 
transformation requires USt < it. To ensure suffi- 
ciently small time discretizations, we meet the condition 
max(C/, e\V\, |eo|i T) < \/6t in all calculations in general. 
After the HS transformation, the remaining fermionic 
Grassmann variables [da, da) appear quadratically and are 
integrated out at the cost of the full path summation over 
the discrete HS Ising spins {s} according to 



{s} o 



(12) 



The full Keldysh GF written in time-discretized (1 < 
k,l < N) form is 

{G-Xi M^V] - {G^aXi+^vS'lr^-^StS^iX^stS^p, 

(13) 
where a, /3 = ± labels the Keldysh branches, and the non- 
interacting GF is 



Goa. 



kl 



duj 

2^ 



e'^'^C^-'^^Goo-H. 



(14) 
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Figure 1 (Color online) Green's function of the non- 
interacting quantum dot in the presence of the leads for 
different settings of gate and bias voltages as well as mag- 
netic field. The inset in (a) shows |ReG'g>r(t — t')| in a 
log-linear representation. 



Note that Go„{t, t') depends only on time differences due 
to time-translational invariance of the noninteracting part 
which holds at thermal equilibrium of the leads. The re- 
spective time discrete version of the self-energy kernel for 
the source term is, cf. Eqs. (|6]l and (|9|, 



kl - 2 



where jp^u 
tm = mSf. 



af3 



af3 
^R,kl 



SmkSa,+ + SmlSl3,+ 



(15) 



lp{tk ~ ti) and the measurement time is 



3 Formulation of the iterative scheme In this sec- 
tion we review the construction for the iterative expressions 
for the Keldysh partition function which form the ISPI 
scheme. For this, we exploit the property (see, e.g., Ref. 
M53ll73ll74l [75l) that each Keldysh component of G^a{t - 
t') in Eq. ([T4]i decays exponentially at long time differ- 
ences {5t\k — l\ — )• oo) for any finite temperature, see Eq. 
( [T4] l. The related time scale is denoted as correlation or 
memory time Tc- In Fig. [T] we show typical examples of 
ReGQ>r(i — t') for different bias voltages V. The expo- 
nential decrease with time is presented in the inset of Fig. 
Sa), where the absolute value |ReGQ^^(t — t')\ is plot- 
d in log-linear representation for the same parameters. 
For large bias voltages and low enough temperatures, e.g., 
at eV > r and T < 0.2r, the decay is superposed by 
an oscillatory behavior Since the lead-induced correlation 
function decays as - cos[eV{t - t')/2]/ sinh[7rT(t - t')], 
the respective correlations decay on a time scale given by 



T^^ ^ max{kBT,eV). Thus, the exponential decay sug- 
gests to neglect lead-induced correlations beyond r^. This 
motivates an iterative scheme which exactly accounts for 
correlations within an interval Tc, while neglecting them 
outside. Notice that the exponential decay is only present 

for finite T and/or V, whereas atT = V = correlations 
die out only algebraically, and our approach is not appUca- 
ble. 



Let us then face the remaining path sum in Eq. ( 12 1. In 
the discrete time representation, we denote by tg = < 
tN — NSt the initial and final times and tk = k5f The 
discretized GF and self-energy kernels for given spin a are 
then represented as matrices of dimension 2N x 2N . For 
explicit calculations, we arrange the matrix elements re- 
lated to Keldysh space (characterized by the Pauli matrices 
t) and to physical times (ife, i;) as r (g) (fc, I). In particu- 
lar, the ordering of the matrix elements from left to right 
(and from top to bottom) follows increasing times. The 
lead-induced correlations thus decrease exponentially with 
growing distance from the diagonal of the matrix. 

For numerical convenience, we evaluate the generating 
function in the equivalent form 

Z[ii]=NY,T{'^^^D,[{s},rf^, (16) 

{s} o 

with D„ — G^^Goa, see Ref. 1531 for details, and thus get 



D 



a0 
M 



;{s},77] = 5apSki+iSt\aGll;uisl 



afi 






<a,kj- 






(17) 

By construction, we have to sum over 2N auxiliary Ising 
spins, which appear line-wise. The total number of possible 
spin configurations is 2^^. 

Next, we exploit the truncation of the GF by setting 
Dki = for \k — l\St > Tc, where Tc = K6t is the corre- 
lation time, with the respective number K of Trotter time 
slices. All GF matrices have a band structure whose band 
width is given by K. Equivalently, we can use for the spin- 
spin correlation the definition 



J 



sf-s^,if\i-j\<K, 
0, else. 



(18) 



We note that in the continuum limit, i.e., K — N,St ^ 
and N — > oo, the approach is formally exact. 

To proceed, we exploit that the determinant of a 

quadratic block matrix £) — is givenby det(D) = 

det(a) det{d — ca ^b). 
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In time space, we obtain the {N x A^)-Keldysh GF band matrix 

(D^^D^^ 



D = ^(i,iv,, 



\ 



£,21 ^22 ^23 



D^^ D"^^ D^^ 



D*^ D*^ 





jjNk-INk 



\ D^Jf^K-i £)NkNk I 

where the single blocks are {K x X)— block matrices defined as (/,?' = 1, ... , Nk) 

D (l-l)K+l,(l' -l)K+l ■ ■ ■ D(i_i)K+l,l'K 



D''' = 



D, 



IK,{1'-1)K+1 



a 



IKI'K 



Without loss of generality, the number N of Trotter slices is chosen as integer such that Nk = N/K and the matrix 
elements D^i follow from Eq. ( |l7| . We keep their dependence on the Ising spins s^ implicit. Each D^i still has a 2 x 
2— Keldysh structure and a 2 x 2spin structure. If we apply the formula for the determinant from above, the generating 
function (fT6| is represented as 



ZM= Y. dct{i?"[s±,...,4]} 






x det {i?(2^^,)[4+i> • ■ • , 4] - ^'M4+i. • ■ - 4k\ [^"[4, • • ■ , 4']] ' ^''[4+1- • • ■ : 4k]] 



(19) 



where the {Nk ~ 1) x [Nk — 1)— matrix D^2.Nk) is obtained from D(i ^Vj^ ) by removing the first line and the first column. 
In order to set up an iterative scheme, we use the following observation: to be consistent with the truncation of the 
correlations after a memory time KSt, we have to neglect terms that directly couple Ising spins at time differences larger 
than Tc, see also Eq. ( [T8] l, consequently matrix products of the form 

(20) 



jji+2,i+i r^i+iJ+11^1 £,(+!,; (j^i.n-''- jji,i+i fjji+ij+i-i-''-jji+ij+2 ^q 



within the Schur complement in each further iteration step. We do not neglect the full Schur complement but only those 
parts which are generated in the second-next iteration step. With this, we rewrite the generating function as 



Nk-1 



ZM^ ^ det{i?"[4,...,4]} n det{i?'+i^'+M4K+i, 



'{l + l)Ki 









1=1 



^ i^{l-l)K+l^- 



,4k] 



D 



i,;+ir„± 



'iif+l' ■ 



'('+ 



i)k]\ 



Exchanging the sum and the product, and reordering the sum over all Ising spins, we obtain 



ZNk [^N-K+1, • ■ • ' ^n\ ' 



(21) 



(22) 






where Zjvk is the last element obtained from the iterative procedure defined by (Z = 1, . . . , Nk — 1) 



Zi+i[s 



± 

IK+l' ■ 



'{l + l)K 



J— 2_^ ^lb'^(l^i)K+l,---,^lK,-^lK+iy-,^(l+l)K]Zn^(l-i)K+l,---,^ 



^k] 



(23) 



(i-l)if+l''"' lA- 



pss header will be provided by the publisher 



The propagating tensor Ai is read off from Eq. (21 1 as 



^,=det{i?'+Wh±,^i,. 
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^^{i+i)k\ 



The iteration starts with Zi [s^ 






(24) 



,4]=det{i?"[s±,...,4]}. 
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Figure 2 (Color online) Extrapolation scheme for raw data 
obtained with ISPl. First the Trotter error, panel (a) and 
then the memory truncation error (b) are extrapolated. No- 
tice that whenever the ISPl data are converged, they are 
numerically exact. 



Figure 3 (Color online) Corrections 51 to the current due 
to finite U as compared Xo U — Q for small to inter- 
mediate on-dot interaction U < F. Other parameters are 
T = O.ir, en = B = 0. We compai-e the ISPl data (red 
symbols) to a 2nd order perturbative calculation, dotted 
Unes guide the eyes. 



The current is numerically obtained by evaluating Eq. 
Q for a small but fixed value of 77; we have taken 77 — 
0.001 for all results shown. The current as a function of 
time I{t,n) shows a transient oscillatory or relaxation be- 
havior at short times. We extract the stationary value / 
when the current has saturated to a plateau. 

3.1 Convergence and extrapolation procedure 

By construction there are two systematic errors in the 
scheme: (i) the Trotter error due to finite time discretiza- 
tion St ~ t/N, and (ii) the memory error due to a finite 
memory time Tc = KSt- The scheme becomes exact in 
the limit K —^ 00 and Jt — >■ 0. We can eliminate both 
errors from the numerical data in the following system- 
atic way: Step 1: We choose a fixed time discretization 
6t and a memory time Tc- A reasonable estimate for Tc 
is the minimum of l/|ey| and 1/T (see above). With 
that, we calculate the current I{St,Tc), and, if desired, the 
differential conductance dI{St,Tc)/dV (the derivative is 
performed numerically for a small AeV = Q.Oir). The 
calculation is then repeated for different choices of St and 
Tc. Step 2: Next, the Trotter error can be eliminated by 
exploiting the fact that it vanishes quadratically for 6t ^ 
[66, 67, 68 1 . For a fixed memory time Tc, we can thus ex- 
trapolate and obtain dI{Tc)/dV — dl{6t — ?> Q,Tc)/dV, 



which still depends on the finite memory time t^. The 
quadratic dependence on St is illustrated in Fig.l2](a) for 
different values of U. Note that each line corresponds to 
the same fixed memory time t^ = 0.5/ F. Step 3: In a last 
step, we eliminate the memory error by extrapolating to 
1/tc — > 0, and obtain the final numerically exact value 
dl/dV = dI{Tc — > oo)/dV . For the dependence on I/tc, 
we empirically find a regular and systematic behavior as 
shown in Fig.|2|b). The Tc —?> 00 value is approached with 
corrections of the order of 1/tc, see Fig.l2tb). 

We have implemented the iterative scheme together 
with the convergence procedure on standard Xeon 2GHz 
machines. Computations are then only possible for K < 7 
due to the limited memory resources available. Typical run- 
ning times for the shown simulation data are approximately 
15 hours for K = 5. 

4 Benchmarking the approach: comparison with 
exact and perturbative results In this section, we dis- 
cuss the results obtained for the Anderson model. We mea- 
sure energies in units of 7^. Unless noted otherwise, all er- 
ror bars for the shown data points, which are due to the 
Trotter and memory extrapolation scheme, are of the order 
of the symbol sizes in the figures. 
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Figure 4 (Color online) Interaction corrections for U = F 
to the nonlinear (main) and linear conductance (inset) cal- 
culated by ISPI (symbols) and by an incoherent rate equa- 
tion as a function of temperature. In the main panel we 
have chosen eV = SF. 



Figure 5 (Color onHne) Comparison of ISPI and RG meth- 
ods in the mixed valence regime. Parameters are chosen as 
U/F = 2 and the temperature is T = OAF. The RG data 
are taken from Ref. Il75l 



We stop short on reporting that we have recovered the 
exact current for the solvable U — case, as a function 
of the bias and the gate voltage. By construction, the ISPI 
method includes all tunneling processes exactly to arbi- 
trary orders. 

In order to benchmark our code for finite U, we com- 
pare the numerical results to a perturbative calculation at 
the charge degeneracy point eq — B — 0, where the in- 
teraction self-energy can be computed up to second order 
in U 1 20 ,"76]. For a detailed comparison, we plot mostly 
the interaction corrections, SA = A{U) — A{U — 0), 
with A being the current /, the linear conductance G, or 
the nonlinear conductance dl/dV, respectively. Figure [s] 
shows the results for SI as a function of the bias voltage 
for U = 0.1F,U = 0.3F and U = F.ForU = O.IF, we 
perfectly recover the perturbative results, which confirms 
the reliability of our code even in the regime of nonlin- 
ear transport. Clearly, the corrections are small and nega- 
tive, which can be rationalized in terms of indications of 
Coulomb blockade physics, as transport is suppressed by 
a finite on-dot interaction. For U = 0.3F, the current de- 
creases even further, and the deviations between the ISPI 
and perturbative results increase. The relative deviation for 
U — 0.3F is already w 30 — 35%, illustrating that pertur- 
bation theory is already of limited accuracy in this regime. 
Although it well reproduces the overall tendency, there is 
a significant quantitative difference. It is even more pro- 
nounced for U = F, as shown in the figure. Here, second- 
order perturbation theory does not even reproduce qualita- 
tive features. 

4.1 Comparison with a master equation ap- 
proach Next, we compare our numerically exact results 
with the outcome of a standard classical rate equation 
calculation II9I I58II . The rate equation is expected to yield 



reliable results in the incoherent (sequential) tunneling 
regime, when F ^ F. Then, a description in terms of 
occupation probabilities for the isolated many-body dot 
states is appropriate. Results for t/ = -T are shown in 
Fig. HI both for the interaction corrections to the nonlin- 
ear and the linear conductance. When the temperature is 
lowered, F < F, quantum coherent interaction effects 
become more important, as seen from the exact ISPI re- 
sults. They are clearly not captured by the master equation 
in the sequential tunneling approximation. However, for 
F ^ F, interaction corrections are washed out, and the 
master equation becomes accurate, cf Fig. |4] Similarly, 
from our ISPI results, we have found (data not shown) that 
interaction corrections are suppressed by an increasing 
bias voltage as well. 

4.2 ISPI vs. tDMRG and fRG The mixed valence 
regime is characterized by the fact that the QD's occupa- 
tion fluctuates between the different charge states, namely 
empty, singly and doubly occupied, see Ref. llTTl for de- 
tails. Furthermore the coupling to the leads is strong, i.e. 
the Kondo temperature F^ ^ F, see below in Sec. 4.3 
When controlling the gate voltage by tuning eo, it is pos 
sible to tune the nanostructure into this limit. From the 
theoretical point of view another energy scale eo enters, 
again, not as a small parameter and simulation methods are 
rare in this regime. For the case of intermediate Coulomb 
repulsion U/F = 2 we have studied |75| the mixed va- 
lence regime and compare our results to other state-of the 
art methods, functional renormalization group (fRG) lfT9]| 
and time dependent density matrix renormalization group 
(tDMRG) 1 2 1,22 1 . The results are shown in Fig.|5]for the 
steady state current I{V) for eg y^ 0, i.e., away from the 
charge degeneracy point. The results from fRG and ISPI 
match perfectly from small bias voltages eV/F w 0.2 
up to the strong non-equilibrium regime. We note that by 
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Figure 6 (Color online) Linear magneto-conductance as 
a function of the gate voltage eo and U = F. The bias 
voltage is chosen as eV = O.OSr and T = O.ir. ISPI 
data (symbols) are connected by dashed lines as guide to 
the eyes. The full lines show the analytical results for van- 
ishing interaction [/ — > 0. Finite magnetic fields split the 
resonant tunneling peak. 



construction, it becomes increasingly cumbersome (and fi- 
nally impossible) to obtain converged ISPI results in the 
limit of vanishing bias voltages and low temperatures f53l, 
since then the correlations do not decay sufficiently to be 
truncated. In the present setup, tDMRG has a tendency to 
overestimate the currents in the mixed valence regime, see 
Ref. ||78l for details, hence we find slight deviations be- 
tween tDMRG and the two other methods (see eV « l.hF, 
lower left panel). Overall agreement away from the sym- 
metric point between the three methods is very good. Fur- 
thermore, I{V) curves show agreement for a wide range of 
magnetic fields, applied to the QD, see Ref. 1751 for details. 

4.3 Small bias regime eV <C -T: For sufficiently 
small bias voltage, the current is linear in V, and we can fo- 
cus on the linear conductance G. Figurel6lshows G'(eo) for 
different magnetic fields B, taking U — F and T — QAF 
(for eV — 0.051^). For B — 0, two spin-degenerate trans- 
port channels contribute, and a single resonant-tunneling 
peak at eo = results. For B y^ 0, the spin-dependent 
channels are split by Aeo — 2B, resulting in a double- 
peak structure. Furthermore, since also U lifts the degen- 
eracy, the spin-resolved levels are now located at eo = 
±{B + U/2) due to the Zeeman splitting. We find an 
interaction-induced broadening, cf. Fig.l6] of the resonant- 
tunneling peak as compared to the noninteracting case. The 
width of the Lorentzian peak profile for i? = is deter- 
mined by F at sufficiently low F, and broadens as F in- 
creases. Here, the double-peak structure, with two clearly 
separated peaks for finite B, is not yet fully developed. The 
two peaks largely overlap, and the distance of the peaks is 
below the expected Aeo — 2B, since tunneling signifi- 
cantly broadens the dot levels. 



Figure 7 (Color online) Log-linear plot of the linear con- 
ductance G as a function of temperature. Above and close 
to the Kondo temperatures Fx, the ISPI simulation results 
(symbols) agree with the results of Ref. 1.79.1 given as full 
lines, see text. 



Next, we address the temperature dependence of the 
linear conductance (numerically evaluated for eV = 
0.05F). In Fig. [t] we show G{F) for different values 
of U (up toU = AF) at eo == B = 0. For U ^ 1.2F, the 
deviations from the U = 0-result is small. For larger U, 
deviations become more pronounced at low temperatures 
where the on-dot interaction becomes increasingly rele- 
vant. Up to present, we have obtained converged results in 
the regime of small bias voltages for interaction strengths 
U < AF for temperatures above or close to the Kondo 
temperature, F > Fk- The corresponding Kondo temper- 
atures are Fk = 0.38r for U ^iF and Fk = 0.293^ for 
U = 4r. In the regime Fk <F < lOT^, we compai-e our 
results to the result of Hamann B77II79I . 



G(T) 



e~ 



1 



HF/Fkh) 



K(r/TKH) + 37r2/4]i/2 



(25) 



for the linear conductance, where Fkh = Fk/1-2, see 
Fig. |7] (solid lines). In Ref. [77J, it has been shown that 
the results of the numerical RG coincide with those of Eq. 
(J25]l in this regime. Fig. IT] illustrates that the agreement 
between the two approaches is satisfactory and shows that 
the ISPI provides reliable results in the linear regime above 
or close to the Kondo temperature. Due to the construction 
of the approach, the situation is more favorable for large 
bias voltages, where short to intermediate memory times 
allow us to obtain convergent results. 

4.4 Large bias regime eV > F: Let us then turn 
to nonequilibrium transport for voltages eV > F. Here, 
the transport window is given by ^ eV, and a double-peak 
structure for dl/dV emerges even for B — 0, see Fig. [S] 
with distance eV between the peaks. We show results for 
eV = 3F but otherwise the same parameters as in Fig. 
[6] For an additional finite magnetic field, each peak of the 
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Figure 8 (Color online) Same as in Fig. p^but for a large 

bias voltage eV = ST. 



2r, eq = B ^ are shown in Fig. l9] Again, as in the 
linear regime, the conductance increases with lower tem- 
peratures, and finally saturates, e.g., at dl/dV ~ e^ /h for 
[/ = and eV = 2F. Clearly the conductance decreases 
when the bias voltage is raised. Increasing U renders this 
suppression yet more pronounced, see also inset of Fig. l9] 
for the corresponding corrections. At high temperatures, 
thermal fluctuations wash out the interaction effects, and 
the interaction corrections die out. 

5 Sustained Franck-Condon blockade in molec- 
ular quantum dots As a second example, we apply 
the ISPI scheme to investigate vibrational effects in the 
nonequilibrium tunneling current through a molecular 
quantum dot, see Ref. ||73]| . We extend the single -impurity 
Anderson model by a linear phonon of frequency i7 (an- 
nihilation operator h) which couples to a single spinless 
electronic level with energy _Eo (operators d/d)). Hence, 
we have the molecular Hamiltonian 
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Figure 9 (Color online) Log-linear representation of the 
temperature dependence of the differential conductance for 
different interaction strengths < U /F < 3 for eV — 
2r, eo = S = 0. 



double-peak structure itself experiences an additional Zee- 
man splitting, resulting in an overall four-peak structure. 
For B = r and the depicted values of U = F, the two 
innermost peaks (closest to eo = 0) overlap so strongly 
that they effectively form a single peak at eo = again. 
The two outermost peaks are due to the combination of the 
finite magnetic field and the bias voltage. 

Increasing the on-dot interaction U downsizes the dif- 
ferential conductance peaks as compared to the nonin- 
teracting case, i.e., the interaction corrections are again 
largest when the level energy matches the chemical poten- 
tial in the leads. Note that the four-peak structure is already 
present in the noninteracting case (with B ^ Q) and, hence, 
is not modified qualitatively by a finite U . 

Finally, we address the temperature dependence of the 
differential conductance dl/dV. ISPI results for eV = 



Hrn = Qb^h+ [Eo + \{h + 6^)] na 



(26) 



with electron-phonon coupling strength A and rici = d^d. 

The short-time propagator on the forward/backward 
branch of the Keldysh contour, e^"'*^. 



^tStH _ 



^^iStHi^^iStiH 



a Trotter breakup, e 

Hi ~ Hjn — flWh, where the auxiliary relation 



then allows for 
^i\ with 



gT«tHi = l-nd 



(27) 
holds. This effectively decouples the electron-phonon in- 
teraction in terms of a three-state variable s^ — 0,±1 
defined at each (discretized) time step tj along the for- 
ward/backward (a — ±) part of the Keldysh contour, 
where 77 = (ij, a). Below, we also use the notation 
77 ± 1 = (^j±i,a) with periodic boundary conditions 
on the Keldysh contour The "Ising spin" variable s^ picks 
up the three terms in Eq. pT] ) and acts like a Hubbard- 
Stratonovich auxiliary field, similar to the Ising field em- 
ployed in the Hirsch-Fye formulation of the Anderson 
model Il53ll66ll67ll68l . The bosonic (phonon) scalar field 
and the fermionic (dot and lead electrons) Grassmann 
fields appearing in the Keldysh path integral are nonin- 
teracting but couple to the time-dependent auxiliary spin 
variable. Hence, those fields can be integrated out ana- 
lytically and the time-dependent current /(im) follows 
from a path summation as above in Sec. l2] The result- 
ing matrix D,,,,' (in time and Keldysh space) depends 
on the complete spin path {s}. Specifically, we obtain 
D — —iB{G^^ — S), where G^^ has spin-dependent 
matrix elements [— iG^^^j . = — s^. We find Sriri' 7^ 
only when s,, = ±1, where it coincides with the usual 
(wide-band limit) expression [80]. Finally, the diagonal 
matrix B (quoted here for eo = 0) with 



B 



vv 



A, 



-X^S^ Y.r,' <^<^'[iG-ph]n,n' + l\'^r,s^i\ 



(28) 
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Figure 10 (Color online) Current / (in units of eF/h) vs 
bias voltage V for the spinless Anderson Holstein model 
for A = O.Sr, ^ = r,Eo = Q, and T = T. The ISPI 
data are depicted as red circles, where the dotted red curve 
is a guide to the eyes only and the error bars are explained 
in the main text. We also show the results of a perturba- 
tion theory in A (solid black curve) and of the rate equation 
(dashed blue curve). The upper (lower) inset shows the cor- 
responding result for T = 3r (T = F/i). 



encapsulates all phonon effects, where Gph is the dis- 
cretized phonon Green's function, see Ref. [81 1, and we 

have used the notation ^0 = landy4±i = ±(l/2)e"^^*t/^. 
We comment shortly on the peculiar convergence prop- 
erties of the present model. Convergence of the extrap- 
olation requires intermediate T or l^-values, for other- 
wise the necessary memory times become exceedingly 
long. For the results below, we have used K < A and 
0.3 < rSt < 0.35. The shown current follows by averag- 
ing over the (5t -window, with error bars indicating the mean 
variance. Additional ISPI runs for 0.18 < r5t < 0.22 and 
0.3 < rSt < 0.4 were consistent with these results, and 
we conclude that small error bars indicate that convergence 
has been reached. When dealing with 3^^'^ summands in 
the evaluation of the generating function in the AH case as 
compared to 2^^ summands for the Anderson model, our 
ISPI code to calculate I {St, K) runs for w 11 CPU hours 
on a 2.93 GHz Xeon processor 

Here, we have once more convinced ourselves that the 
numerical ISPI results for the I{V) curves are consistent 
with known analytical theory in the respective parameter 
limits. We employ the following perturbative methods: (i) 
For X/r ^ 1, perturbation theory in the electron-phonon 
coupling applies and yields a closed KV) expression for 
arbitrary values of all other parameters f82l . We note that 
the solution of the AH model with a very broad dot level 
1 83 84 1 corresponds to this small- A regime, (ii) For high 
temperatures, T ^ r,a description in terms of a rate equa- 
tion is possible |80|. We here use the sequential tunneling 
approximation with golden rule rates II851 . For small A, the 



Figure 11 (Color online) Same as Fig. 10 but for Q = 
O.bF and X = F. The main panel is for F = F and com- 
pares the ISPI results to NEBO predictions. The insets are 
for T = 3F and T = F/3, respectively, where also the rate 
equation results are shown. Notice that in contrast to ISPI, 
the rate equation predicts an unphysical current blockade 
for T = F/3. 



corresponding results match those of perturbation theory, 
while in the opposite strong-coupling limit, the Franck- 
Condon blockade occurs and implies a drastic current sup- 
pression at low bias voltage |30,86|. (iii) For small os- 
cillator frequency, i7 <C min(I^, eV), the nonequilibrium 
Born-Oppenheimer (NEBO) approximation is appropriate 
and allows us to obtain I{V) from a Langevin equation 
for the oscillator r87'"881. For small A, this approach is 
also consistent with perturbative theory, while for high F, 
NEBO and rate equation results are found to agree. For 
clarity, we focus on a resonant level with E^ = here. 
The case of weak electron-phonon coupling, A ~ Q.5F is 
shown Fig. [To] We compare our ISPI data f or i? := I^ to the 
results of perturbation theory in A and of the rate equation. 
Perturbation theory essentially reproduces the ISPI data. 
The rate equation is quite accurate for high temperatures, 
but quantitative agreement with ISPI was obtained only for 
F > lOF. We note that the ISPI error bars increase when 
lowering F due to the growing memory time {Tc) demands. 



Next, Fig. 11 shows ISPI results for a slow phonon 
mode, n = rj2, with larger electron-phonon coupling 
A = -T. In that case, perturbation theory in A is not reli- 
able and likewise, the rate equation is only accurate at the 
highest temperature (T = 3F) studied, cf. the upper left 



inset of Fig. 1 1 However, we observe from Fig. 1 1 that for 



such a slow phonon mode, NEBO provides a good approx- 
imation for all temperatures and/or voltages of interest. We 
conclude that the ISPI technique is capable of accurately 
describing three different analytically tractable parameter 
regimes. 

In the limit of strong electron-phonon coupUng A, the 
classical rate equation predicts a Franck-Condon block- 



Copyright line will be provided by the publisher 



12 



S. Weiss et al.: Iterative summation of path integrals 






' 


1 


/ 


' 


1 




■ T 


= r 


/ 






■ 


- 


/ 








; 
















/ 




^ 




. .- 


/ 




^ 


^ 


■ ^ 


- 














<>■•'' 


r-T 


- 


"] 







2 


4 




6 


> 



A = o.5r T=o.2r 
x=i.5r 
>. = 2.5r 




Figure 12 (Color online) ISPI data for the liy) curves 
for the spinless Anderson Holstein model from weak (A — 
O.S-T) to strong (A = 41^) electron-phonon coupling, with 
Q = 2r. The main panel is for T — 0.21^, the inset for 
T = r. We used a dense voltage grid yielding smooth 
I{V) curves. Error bars are not shown but remain small, 
cp. Fig. [To 



ade of the current for low bias and T ^ F 1.86 J . Suf- 
ficiently large A can be realized experimentally, and the 
Franck-Condon blockade has indeed been observed in sus- 
pended carbon nanotube quantum dots |30|. For a nonequi- 
librated phonon with intermediate-to-large A, understand- 
ing the Franck-Condon blockade in the quantum coherent 
regime of low temperature, T < I^, is an open theoretical 
problem. Here, multiple phonon excitation and deexcita- 
tion effects generate a complicated (unknown) nonequilib- 
rium phonon distribution function, and the one-step tun- 
neling interpretation in terms of Franck-Condon factors 
between shifted oscillator parabolas |86| is no longer ap- 
plicable. We here study this question using ISPI simula- 
tions, which automatically take into account quantum co- 
herence effects. In Fig. 12 the crossover from weak to 



strong electron-phonon coupling A is considered. The inset 
shows I{V) curves for T — F, where we observe a current 
blockade for low voltages once A > 2F. The blockade be- 
comes more pronounced for increasing A and is lifted for 
voltages above the polaron energy A^ / i7 1 86 1 . Remarkably, 
the Franck-Condon blockade persists and becomes even 
sharper as one enters the quantum-coherent regime (here, 
T = Q.2F), despite of the breakdown of the sequential tun- 
neling picture. We also observe a nonequilibrium smearing 
of phonon step-like features in the I{V) curves in Fig. 12 
cf. alsoRefs. EOllMl- 

6 Nonequilibrium quantum dynamics in the mag- 
netic Anderson model The third example to which we 
have applied the ISPI scheme |74| is the magnetic Ander- 
son model. We extend the Hamiltonian Eq. ([T]) in the pres- 
ence of Coulomb interactions, see Eq. ( [TO] i, by a magnetic 



impurity localized on the quantum dot which interacts via 
an exchange interaction with the spins of the confined elec- 
trons on the QD. The magnetic part of the Hamiltonian 
reads 



H, 



imp 



^int — 



Ai. 



+ jT^{d\.d^ - d\d^) + -{T+d\d^ + T^d\.d^) . 



HI 



Hi 



(29) 



The generating function Z[rj\ is obtained by integrating 
again over the corresponding Grassmann fields for dot and 
lead operators as well as the discrete paths {r} and {(,} for 
the real spin variables and the HS Ising fields, respectively, 
and 



ZM 



E 



T^[CkpadaCkpada]{-lf 



iJSt 
2 



P[{T}]e 

(30) 

The path sums over impurity and HS spin-fields are per- 
formed over the 2A^-tuples {tj} = {t2n,---,ti) and 
{Q} = {C2N, ■ • • , Ci) with Tj,(j = ±1. Within an impu- 
rity path {r}, m flip-flop transitions occur on the Keldysh 
contour, where i of them lie on the lower branch. The ac- 
tion S includes tunneling and lead effects as in Eq. (|6|. 
Correspondingly, the magnetic part of the action is 



iS 



S, 



^imp^t 



N 



imp 



fe=2 



{Tk — T2N-k+ 



(31) 



The polynomial ^^[{t}] in Eq. ( pOl l depends on the impu- 
rity path {r} = {■'"~''}({t~}) for the forward (backward) 
branch of the contour Then, we collect all indices of 
the flips into the tuple F^ = (fc^_^, . . . , fcj^) (sorted 
in ascending order) along the forward path {r+} := 
(rjv, . . . , Ti) with Tfe+ ^ Tfc+_i for afl k+ £ F+^. Accord- 
ingly, TjT = {kj , . . . ,k^)isthe tuple of ascending flip in- 
dices along the backward path {r~} := {t2n, ■ ■ ■ , tn+i) 
with Tj,- y^ Tf.-+i for all fc^ e T^jj" . Note that a flip index 
on the backward path is labelled according to the smaller 
step index of the flipping spins corresponding to the later 
time. The impurity polynomial can be expressed in terms 
of the electronic Grassmann fields as 



p[{r}] := n ^'dir, n ^' 



rjr. 



■fc-1 



(32) 



jeT„. 



feeT„t 



Figure [13] illustrates an example of an impurity path. Col- 
lecting all pieces, the remaining formally exact expression 
for the Keldysh generating function is 

2[V]= 5](P[{r}])ndet{(*Gf[{T,C},r,])-i}. (33) 

{rX} 
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Tit) 



1+ ,2^ ,3+ ,4+ ,5+ ,6^ ,7+ ,8^ ,9+ 




14 13 12 11 10 



Figure 13 Exemplary impurity Keldysh path (blue line). 
The Keldysh contour is divided into iV — 1 = 8 segments 
of length St between 2N — 18 time vertices. The impurity 
path (tuple of black and red arrows) realizes m = A flip- 
flops along the contour 



The Keldysh partition function is given as a sum over ex- 
pectation values of the polynomial P of Grassmann num- 
bers in a system with Green's function G^*. In passing, we 
note that it is possible to express the expectation values for 
the polynomials P in terms of Green's functions for the 
interacting system. The details are given in Ref. 1(741 . Af- 
ter applying Wick's theorem we obtain explicit expressions 
for the polynomial |74J. Using matrix elements {Scr)k.i — 
—i{d%''<rj) — {C^)q^rn the final expression for the gen- 
erating function follows as 



1 

0.8 
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Figure 14 (a) The expectation value (r^) of the impurity 
orientation as a function of propagation time for different 
strengths J of an anti-ferromagnetic electron-impurity in- 
teraction. The initial preparation of the system at t = — oo 
is spin-up {tz{Q) = ti — 1], for the other parameters see 
text. The polarization decays faster with increasing J. (b) 
Impurity relaxation rate r^^ as a function of the exchange 
coupling strength. Again the relaxation is faster with in- 
creasing coupling strength. In both panels the temperature 

isT = r. 



Rapidly decreasing weights of the paths may not be 
(over-)compensated by their increasing numbers for < 
rrij < K — 1, since each contribution is small and the 
number of paths decreases again for larger nij > K. 



xp{^^imp} IT dct z(Gf )-i det H„, (34) The behavior of the impurity weights is illusti-ated as fol- 



X c: 



where the summation over impurity paths is restricted to 
tuples {t} with ti = T2n ~ ^i, i.e., correct boundary con- 
ditions along the Keldysh contour are fulfilled. The limit 
(5f — > appears explicitly here, since there is no continu- 
ous measure used for the discrete spin paths, neither for the 
HS- nor for the impurity spins. In order to reduce the ex- 
ponentially growing number of contributing paths (^ 4^^, 
due to real spin and HS spins) without affecting the accu- 
racy, we may exploit that the propagating tensor depends 
on the number rrij of flip-flops in path segment j and 
< rrij < 2{K — 1) along the Keldysh contour We ob- 
serve that the weight of each segment is smaller, the more 
flip-flops it contains. On the other hand, the number of path 



2(^-1) 



with 



segments {t}j with rrij flip-flops (given by 4(7, 
CJ^ — nl/[kl{n — k)l]) grows as long as < nij < K — I, 
but decreases again when K < rrij < 2{K — 1). As a con- 
sequence, for any observable there exists a maximal m™"^ 
such that contributions from paths with rrij > m™"^ < 
2{K — 1) could safely be disregarded in the numerical iter- 
ation. Of course m™"^ is chosen depending on the model 
parameters and the observable under investigation. 



lows. Consider the case when rrij is close to the max- 
imum 2{K — 1). Both path classes with nij — and 
rrij — 2{K — 1) contain the same number of elements 
(four), while each path contribution in the second class is 
weighted by {J5t/2y^^^^\ For typical values of K — 4, 
StP = 1/2, and J ^ F, the weight is - 2.5 x lO^*. This 
also holds for afl K < rrij < 2{K - 1). Since to™" is 
unknown a priori, we include it into our code as an addi- 
tional parameter Then, we perform a numerical estimate 
by a spot sample of the parameter space. It turns out that 
for the considered cases, it is sufficient already to choose 
TO™" = 2. This drastically reduces the CPU running times 
from more than one month to typically three to five days. 

6.1 Impurity dynamics We focus on tiansport fea- 
tures caused by the magnetic impurity in this section. We 
emphasize that novel dynamical and transport features are 
mediated by the transverse or flip-flop interaction iJj^j, 
given in Eq. ( |29l ). Without the possibility of flip-flops the 
orientation of the impurity spin and its quantum state could 
not change. The remaining longitudinal part Hr^^ of the 
interaction causes a renormalization of rates and energies 
which appears as effective magnetic field. Necessarily, flip- 
flop processes are involved from the beginning to investi- 
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gate the non-trivial impurity dynamics by considering the 
time dependence of the impurity orientation (r^). In all 
presented results below, the impurity is fully polarized at 
t ^ —oo and the coupling to the leads is switched on. 



In Fig. 14 (a) we present the time evolution of the 
impurity polarization (r^) for different values of the ex- 
change interaction J. The remaining parameters are <l>jj = 
A = Z\in,p = [/ = 0, and T = r and eV = 0.6^. 
The impurity polarization shows a clear exponential de- 
cay {Tz)(t) oc e"*-*"*'-*"^!* , well described by a single re- 
laxation rate for intermediate to long propagation times. A 
faster decay is observed as the impurity interacts stronger 
with the electron spins. The parameters are chosen to yield 
an isotropic (symmetric with respect to [relative] spin ori- 
entations) model system. In this case the antiferromagnetic 
interaction favors antiparallel orientation of electron- and 
impurity spin. Over long propagation times, the coupling 
to the unpolarized leads then destroys any polarization of 
the impurity. It is therefore reasonable to assume, that the 
rates for up- and down flips are equal. While the impurity 
interaction energy is comparable to the tunneling coupling 
and considerably affects the transport behavior as we show 



below (see Fig. 15 i, the rather high temperature and bias 



voltage nevertheless reduce the relevance of coherent dy 
namics due to on-dot interactions to a secondary role 

In Fig. 



14 



:b) 



We next investigate the relaxation rate t^ 
for T = r. We present results for varying J and U — 
and three different bias voltages. These show a nearly 
quadratic behavior growing from zero (no relaxation) in 
the sense that for a fit of the results for < J < r/2 
to a polynomial function aJ^ the exponent b lies between 
^ 1.8 and ^ 1.9. An exact quadratic dependence of t^^ 
on J is obtained only when the dynamics is strongly dom- 
inated by sequential (incoherent) flip-flop processes 1741 . 
This is only realized when J -^ F. A sequential flip- 
flop process consists of three elementary components: 
the actual flip-flop and two tunnelling processes of single 
electrons with opposite spins (not necessarily in that or- 
der). Since they evolve coherently, these components form 
an effective spin-flip process \x,t) — > \x,—t), where 
X G {0, a, d} and the underlying flip-flop nature is masked 
by the tunnelling electrons. 

6.2 Charge current for finite impurity interaction 
and Coulomb repulsion In the deep quantum regime, 
where no small parameter exists, ISPI is able to describe 
physical properties not predictable by perturbative meth- 
ods. In this section, we study how the current behaves as 
functions of bias voltage. Coulomb interaction and tem- 
perature, respectively. 

Fig. 15 (a) shows that the flip-flop term iJj^j has a con- 
siderably smaller influence on the charge current atT = F 
(incoherent regime) than the longitudinal part of the in- 
teraction in Eq. (|29]l. Despite the qualitatively similar be- 
havior of the Landauer-Biittiker (LB) current and the exact 
data, the flip-flop scattering causes an additional significant 
current drop that grows for growing J. A finite Coulomb 
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Figure 15 (a) Charge current as a function of the ex- 
change interaction J for two values of the on-dot interac- 
tion t/ at r = r and eF = 0.6/". The solid line marks the 
Landauer-Biittiker result, where exchange correlations are 
treated on a mean field level, see Ref 1741 . (b) Comparison 
of (i) the LB current (solid lines), (ii) the Coulomb inter- 
acting current without flip-flop scattering ("no flips", red 
circles), (iii) the current without Coulomb scattering but 
full impurity interaction ("mean-field C/", green diamonds, 
see the text for explanation), and (iv) the fully interacting 
current ("full int.", blue squares) in their dependence on the 
Coulomb interaction U for T = F . The other (non-zero) 
parameters are J = I" and eV = 2F. 



interaction of U = F/2 increases the resistivity of the dot 
and the ISPI points are consistently lower than the LB val- 
ues. The voltage is chosen as eV — 0.6F. In Fig 
foTT = F four different current curves are shown 



151 (b), 
one 

for each possibility to either have (i) only mean field dy- 
namics, regarding J (LB), (ii) the full Coulomb interac- 
tion without flip-flop processes ("no flips"), (iii) flip-flop 
dynamics without Coulomb fluctuations ("mean-field IF'), 
and (iv) the fully interacting dot ("full int."). For J — F 
and V — 2F, the Coulomb energy is varied between 
< U < F. The situation "mean-field [/" is implemented 
by setting ^^ = U/2 and the HS parameter A = to il- 
lustrate the effect of the "classical" part of the Coulomb 
interaction. Only for the "single-interaction" currents ("no 
flips"), we show the error bars. We do not show a margin 
of confidence for the fully interacting case in order that the 
error data remain comparable. Calculating the "full int." 
current is a time consuming task and thus, the extrapola- 
tion involves considerably fewer data points. Nevertheless, 
this does not render these values unreliable (we still see a 
compelling linear behavior of the l/r^ extrapolation with 
errors of the order of 1% based on the sample standard de- 
viation). Both the mean-field current and the current with- 
out Coulomb scattering show only a weak dependence on 
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U due to the single-particle energy shift. The current with 
full Coulomb interaction but fixed impurity shows a local 
maximum for U ^ r/2. In this case, the fixed impurity 
acts as an effective static magnetic field. The ISPI values 
for the fully interacting dot vary strongly over the consid- 
ered U interval, but are scattered around the "no flips" and 
"mean-field U" curves. 

As long as the Coulomb interaction is small, differ- 
ences in the data arise from including or excluding flip-flop 
processes. Hence, the rather good agreement of the U = 
values suggests that even at this, temperature flip-flop pro- 
cesses alone affect the current only weakly. Nevertheless, 
for decreasing temperatures, the flip-flop processes start to 
influence the current more strongly, which results in an in- 
creased resistivity. The case of the "no flip" current (fixed 
impurity) is equivalent to a Coulomb-interacting single- 
level quantum dot in a magnetic field. This effect is caused 
by the broadening of the dot's joint density of states due to 
the Coulomb fluctuations. 

7 Conclusions In summary, we have reviewed a 
scheme for the iterative summation of real-time path in- 
tegrals (ISPI) and applied it to prototypical problems of 
quantum transport through an interacting quantum dot 
coupled to metallic leads held at different chemical poten- 
tials. After integrating over the leads' degrees of freedom, 
a time-nonlocal Keldysh self-energy arises. Exploiting the 
exponential decay of the time correlations at finite tem- 
perature allows us to introduce a memory time Tc beyond 
which the correlations can be truncated. Within Tc, corre- 
lations are fully taken into account in the corresponding 
path integral for the Keldysh generating function. Then, 
through a discrete Hubbard-Stratonovich transformation, 
interactions are transferred to an auxiliary spin field, and 
an iterative summation scheme is constructed. The remain- 
ing systematic errors due to the finite time discretization 
and the finite memory time Tc are eliminated by a refined 
Hirsch-Fye-type extrapolation scheme, rendering the ISPI 
numerically exact. 

The scheme has been applied to the canonical exam- 
ple of a single-impurity Anderson dot with Coulomb in- 
teraction U. This allows us to carefully and systematically 
check the algorithm. For linear transport, we have recov- 
ered results from second-order perturbation theory in U in 
the limit of very small interaction strength, but found sig- 
nificant deviations already for small-to-intermediate values 
of U. In the incoherent sequential regime, we recover re- 
sults from a master equation approach. We have further- 
more reproduced the linear conductance above the Kondo 
temperature. In addition, we have investigated the regime 
of correlated nonlinear transport, where, in our opinion, 
the presented method is most valuable. The nonequilibrium 
Kondo regime, representing an intermediate-to-weak cou- 
pling situation, seems tractable by the ISPI scheme. 

Our approach is, in fact, similar in spirit to the well- 
established concept of the quasi-adiabatic path integral 



(QUAPI) scheme, introduced by Makri and Makarov 118911 
in its iterative version. This method has been developed 
to describe the dynamics of a quantum system coupled to 
bosonic environments, see also Refs. II70II90I . 

In a second example, we have applied the ISPI tech- 
nique to the spinless Anderson-Holstein model as well, 
which is the simplest nonequilibrium model for molecu- 
lar quantum dots with a phonon mode. Our formulation 
exploits a mapping to an effective three-state system and 
reproduces three analytical theories valid in different pa- 
rameter regions. This extension of the ISPI approach then 
captures the full crossover between those limits. For strong 
electron-phonon coupling and a nonequilibrated phonon 
mode, we find that the Franck-Condon blockade becomes 
even more pronounced as one enters the deep quantum co- 
herent regime. 

The complex system of an incorporated spin- 1/2 mag- 
netic impurity in a quantum dot has been subject of a third 
study, focusing on the real-time dynamics in the presence 
of Coulomb interactions. We include the impurity inter- 
action on the same level as the other interactions, which 
results in an additional sum over impurity paths. An ef- 
ficient truncation scheme nevertheless provides accurate 
results for the coupled spin dynamics. Results are given 
for a quantum spin- 1/2 impurity on the dot, whereas the 
generalization to an impurity with a larger spin is possi- 
ble. For a small impurity interaction, where sequential flip- 
flops dominate the impurity dynamics, we have found good 
agreement with a classical rate equation, see Ref |74| as 
well. This is a useful tool to gain insight into the dom- 
inating processes in the incoherent regime. Relaxation is 
described reasonably well by a rate equation when lead- 
induced coherences are absent. 

In the deep quantum regime, however, we find that the 
ISPI method is the only tool to obtain both the correct order 
of magnitude and the qualitative features of the relaxation 
rate as it depends on the system parameters U and J. The 
same holds for the influence of J and U on the current in 
this interesting corner of the parameter space. Most impor- 
tantly, the ISPI scheme proves to be useful to cover the 
full cross-over regime where no small parameter exists and 
thus any perturbative approach becomes invalid. 

We have provided a first glimpse on the interesting 
new physics that comes into reach with the ISPI scheme. 
Compared to other approaches, it has several advantages 
(e.g., numerical exactness, direct nonequilibrium formula- 
tion, no sign problem), but is, on the other hand, computa- 
tionally more costly than most other techniques, especially 
for strong correlations and/or low energy scales (tempera- 
ture, voltage). 
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